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Models: blow-up and steady states 
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Abstract 



Nonlinear Noisy Leaky Integrate and Fire (NNLIF) models for neurons networks can be written 
as Fokker-Planck-Kolmogorov equations on the probability density of neurons, the main parameters 
in the model being the connectivity of the network and the noise. We analyse several aspects of 
the NNLIF model: the number of steady states, a priori estimates, blow-up issues and convergence 
toward equilibrium in the linear case. In particular, for excitatory networks, blow-up always occurs 
for initial data concentrated close to the firing potential. These results show how critical is the 
balance between noise and excitatory /inhibitory interactions to the connectivity parameter. 



Key-words: Leaky integrate and fire models, noise, blow-up, relaxation to steady state, neural 
qq ' networks. 

d 



AMS Class. No: 35K60, 82C31, 92B20 



1 Introduction 

The classical description of the dynamics of a large set of neurons is based on deterministic/stochastic 
^ differential systems for the excitatory-inhibitory neuron network \12\ [23] . One of the most classical 

models is the so-called noisy leaky integrate and fire (NLIF) model. Here, the dynamical behavior of 
the ensemble of neurons is encoded in a stochastic differential equation for the evolution in time of the 
averaged action potential of the membrane v(t) of a typical neuron representative of the network. The 
neurons relax towards their resting potential Vi in the absence of any interaction. All the interactions 
of the neuron with the network are modelled by an incoming synaptic current I{t). More precisely, 
the evolution of the action potential follows, see [U [U [20j [7] 

C m ^ = -9L(V-V L )+I(t) (1.1) 

where C m is the capacitance of the membrane and gL is the leak conductance, normally taken to be 
constants with r m = gi/Cm ~ 2ms being the typical relaxation time of the potential towards the leak 
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reversal (resting) potential Vl ~ —70mV. Here, the synaptic current takes the form of a stochastic 
process given by: 

Ce Ci 

j=l j i=l j 

where <5 is the Dirac Delta at 0. Here, Je and J/ are the strength of the synapses, Ce and Cj are 
the total number of presynaptic neurons and t l E j and t\- are the the times of the j^-spike coming 
from the i^-presynaptic neuron for excitatory and inhibitory neurons respectively. The stochastic 
character is embedded in the distribution of the spike times of neurons. Actually, each neuron is 
assumed to spike according to a stationary Poisson process with constant probability of emitting a 
spike per unit time v. Moreover, all these processes are assumed to be independent between neurons. 
With these assumptions the average value of the current and its variance are given by fj,c = bv with 
b = CeJe — CjJj and a c = (CeJe + CjJj)u. We will say that the network is average-excitatory 
(average-inhibitory resp.) if b > (b < resp.). 

Being the discrete Poisson processes still very difficult to analyze, many authors in the literature 
[U HI [2"m [T3] have adopted the diffusion approximation where the synaptic current is approximated by 
a continuous in time stochastic process of Ornstein-Uhlenbeck type with the same mean and variance 
as the Poissonian spike-train process. More precisely, we approximate I(t) in (|1.2p as 

I(t) dt « fj, c dt + ac dB t 

where Bt is the standard Brownian motion. We refer to the work [20] for a nice review and discussion 
of the diffusion approximation which becomes exact in the infinitely large network limit, if the synaptic 
efficacies Je and Ji are scaled appropriately with the network sizes Ce and Cj. 

Finally, another important ingredient in the modelling comes from the fact that neurons only fire 
when their voltage reaches certain threshold value called the threshold or firing voltage Vp ~ —50mV. 
Once this voltage is attained, they discharge themselves, sending a spike signal over the network. We 
assume that they instantaneously relax toward a reset value of the voltage Vr ~ — 60mV. This is 
fundamental for the interactions with the network that may help increase their action potential up 
to the maximum level (excitatory synapses), or decrease it for inhibitory synapses. Choosing our 
voltage and time units in such a way that C m = ql = 1, we can summarize our approximation to the 
stochastic differential equation model (jl.ip as the evolution given by 

dV = {-V + V L + n c )dt + a c dB t (1.3) 

for V < Vf with the jump process: V(t„ ) = Vr whenever at to the voltage achieves the threshold 
value V(t~) = Vf; with Vl < Vr < Vf- Finally, we have to specify the probability of firing per unit 
time of the Poissonian spike train v. This is the so-called firing rate and it should be self-consistently 
computed from a fully coupled network together with some external stimuli. Therefore, the firing rate 
is computed as v = v ex t + N(t) where N(t) is the mean firing rate of the network. The value of N(t) 
is then computed as the flux of neurons across the threshold or firing voltage Vp- We finally refer to 
|llj for a nice brief introduction to this subject. 

Coming back to the diffusion approximation in (|1 .3|) , we can write a partial differential equation for 
the evolution of the probability density p(v,t) > of finding neurons at a voltage v E (— oo, Vf] at a 
time t > 0. Standard Ito's rule gives the backward Kolmogorov or Fokker-Planck equation 

%(v, t) + £ [h(v, N(t))p(v, t)] - a [N(t))Pz(v, t) = 5(v - V R )N(t), v<V F , (1.4) 
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with h(v,N(t)) = —v + Vl + fJ> c an d a(N) = We have the presence of a source term in the 

right-hand side due to all neurons that at time t > fired, sent the signal on the network and then, 
their voltage was immediately reset to voltage Vr. Moreover, no neuron should have the firing voltage 
due to the instantaneous discharge of the neurons to reset value Vr, then we complement (jl.4p with 
Dirichlet and initial boundary conditions 

p(V F ,t) = 0, p(-oo,t) = 0, p(v,0)=p°(v). (1.5) 

Equation (II, 4ft should be the evolution of a probability density, therefore 

V F rVp 

p(v, t)dv = I p°(v) dv = 1 



for all t > 0. Formally, this conservation should come from integrating (jl,4p and using the boundary 
conditions (|1.5p . It is straightforward to check that this conservation for smooth solutions is equivalent 
to characterize the mean firing rate for the network N(t) as the flux of neurons at the firing rate voltage. 
More precisely, the mean firing rate N(t) is implicitly given by 

N(t):=-a{N(t))^(V F ,t)>0. (1.6) 

Here, the right-hand side is nonnegative since p > over the interval [— oo, Vf] and thus, ^{Vf, t) < 0. 
In particular this imposes a limitation on the growth of the function N i— > a(N) such that (jl.6p has a 
unique solution N. 

The above Fokker-Planck equation has been widely used in neurosciences. Often the authors prefer 
to write it in an equivalent but less singular form. To avoid the Dirac delta in the right hand side, 
one can also set the same equation on (— oo, Vr) U (Vr, Vf] and introduce the jump condition 

p(V R -,t)=p(V+,t), ^p(VR,t) - §f0#,t) = N(t). 

This is completely transparent in our analysis which relates on a weak form that applies to both 
settings. 

Finally, let us choose a new voltage variable by translating it with the factor Vl + bv ex t while, for 
the sake of clarity, keeping the notation for the rest of values of the potentials involved Vr < Vf- In 
these new variables, the drift and diffusion coefficients are of the form 

h(v,N) = -v + bN, a(N) = a + ai N (1.7) 

where b > for excitatory-average networks and b < for inhibitory-average networks, ao > and 
0,1 > 0. Some results in this work can be obtained for some more general drift and diffusion coefficients. 
The precise assumptions will be specified on each result. Periodic solutions have been numerically 
reported and analysed in the case of the Fokker-Planck equation for uncoupled neurons in |16l \17\ . 
Also, they study the stationary solutions for fully coupled networks obtaining and solving numerically 
the implicit relation that the firing rate iV has to satisfy, see Section 3 for more details. 

There are several other routes towards modeling of spiking neurons that are related to ours and 
that have been used in neurosciences, see [9]. Among them are the deterministic I&F models with 
adaptation which are known for fitting well experimental data [3]. In this case it is known that 
in the quadratic (or merely superlinear) case, the model can blow-up [22]. One can also introduce 



3 



gating variables in neuron networks and this leads to a kinetic equation, see [6] and the references 
therein. Another method consists in coding the information in the distribution of time elapsed between 
discharges [HHIK!, this leads to nonlinear models that exhibit naturally periodic activity, but blow-up 
has not been reported. 

In this work we will analyse certain properties of the solutions to (|1.4p -( fL~5j) with the nonlinear term 
due to the coupling of the mean firing rate given by (jl.6p . Next section is devoted to a finite time 
blow-up of weak solutions for (11.4p - (11.6p , In short, we show that whenever the value of b > is, we can 
find suitable initial data concentrated enough at the firing rate such that the defined weak solutions 
do not exist for all times. This implies that this model encodes complicated dynamics. As long as 
the solution exists in the sense specified in Section [21 we can get apriori estimates on the L^-norm 
of the firing rate. Section [3] deals with the stationary states of (|1.4p — (|1.6j) . We can show that there 
are unique stationary states for b < and a constant but for b > different cases may happen: one, 
two or no stationary states depending on how large b is. In Section [4j we discuss the linear problem 
6 = with a constant for which the general relative entropy principle applies implying the exponential 
convergence towards equilibrium. Finally, Section [5] is devoted to some numerical simulations of the 
model showing some of the results here and getting some conjectures about the nonlinear stability of 
the found stationary states. 



2 Finite time blow-up and apriori estimates for weak solutions 



Since we study a nonlinear version of the backward Kolmogorov or Fokker-Planck equation (jl.4p , we 
start with the notion of solution: 



Definition 2.1 We say that a pair of nonnegative functions (p,N) with p £ L c 



-oo 



N G Li 



,Vf)), 



loc,- 



is a weak solution of (jl.4f) (jl.7() if for any test function cp(v, t) G C°°((— oo, Vf] x [0, T]) 



such that G L°°((-oo,V F ) x (0,T)), we have 



T r V, 



p(v,t) 



3 $ d< t>Ut AT\ 9 " ^ 

~dt ~ ~~dv ~ a ~du* 



dv dt 



N{t)[<j>{V R ,t) - <j>(V F ,t)]dt 



(2.1) 



+ 



Vf 



p°(v)<j){0,v)dv 



V F 



p(v, T)(f>(T, v) dv. 



Let us remark that the growth condition on the test function together with the assumption (|l,7p 
imply that the term involving h(v,N) makes sense. By choosing test functions of the form i[)(t)(f>(v), 
this formulation is equivalent to say that for all (j)(v) G C°°((— oo, Vp]) such that v^ G L°°((— oo, Vf)), 
we have that 



d 



— J (p(v)p(v,t)dv 



p(v, t)dv + N(t)[cj>(V R ,t) - <j>(V F ,t)] (2.2) 



holds in the distributional sense. It is trivial to check that weak solutions conserve the mass of the 
initial data by choosing <fi = 1 in (|2.2p . and thus, 



/Vf [Vf 
p(v,t)dv= / p°(v)dv 
-oo J — oo 



1 . 



(2.3) 



The first result we show is that global-in-time weak solutions of (jl.4j) - (jl.6j) do not exist for all initial 
data in the case of an average-excitatory network. This result holds with less stringent hypotheses on 
the coefficients than in (|1.7p with an analogous notion of weak solution as in Definition 12.11 
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Theorem 2.2 (Blow-up) Assume that the drift and diffusion coefficients satisfy 

h(v, N) + v>bN and a(N) > a m > 0, (2.4) 

for all —oo < v <Vf and all N > 0, and let us consider the average- excitatory network where b > 0. 
// the initial data is concentrated enough around v = Vf, in the sense that 

f v F 

/ e^W) dv 

J — oo 

is large enough with fi > max(^, ^), then there are no global-in-time weak solutions to ()1.4|) fjl.6[) . 

Proof. We choose a multiplier 4>(v) = e 1 " 3 with /i > and define the number 

A = <f>{V F ) - 4>(V R ) > Q 
bfj, 

by hypotheses. For a weak solution according to (|2.ip . we find from (|2.2p that 



df 



d r-V F rV F rV F 

4>(v)p(v,t)dv > ill (bN(t) - v)<j)(v)p(v,t)dv + /u 2 a m / <j>(v)p(v,t) dv - \bfxN(t) 

X) J -co J —oo 

rV F 

>fj, 4>(v)p(v, t) dv [bN(t) + /ua m - V F ] - XfibN(t) (2.5) 

J — oo 

where (|2.4h and the fact that v S (— oo,Vp) was used. Let us now choose large enough such that 
fj,a m — Vf > according to our hypotheses and denote 

fV F 

M„(t) = / <f>(v)p(v,t)dv, 

J —oo 

which satisfies 

j t M»(t) >bnN{t)[M^t)-\]. 

If initially M^(0) > A and using Gronwall's Lemma since N(t) > 0, we have that M^{t) > A, for all 
t > 0, and back to (|2.5p we find 



dt 

which in turn implies, 



4>(v)p(v, t) dv > n(fj>a m - Vf) / 4>(v)p(v, t) dv 



v F rV F 

(f>{v)p{v, t) dv > e^ am - VF)t / <f>(v)p°(y) dv. 



On the other hand, since p(v,t) preserves the mass, see (12.30 . and fi > then 



{v)p(v,t)dv < e^ Vp , 



leading to a contradiction. 
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It remains to show that the set of initial data satisfying the size condition in the statement is not 
empty. To verify this, we can approximate as much as we want by smooth initial probability densities 
an initial Dirac mass at Vf which gives the condition 



e» VF > A 



together with fia m > Vf- 



This can be equivalently written as 



1 _ e -Kv F ~v R ) 



and fi > — . 



b 



Choosing /i large enough, these conditions are obviously fulfilled. □ 

As usual for this type of blow-up result similar in spirit to the classical Keller-Segel model for 
chemotaxis [21 [8], the proof only ensures that solutions for those initial data do not exist beyond a 
finite maximal time of existence. It does not characterize the nature of the first singularity which 
occurs. It implies that either the decay at infinity is false, although not probable, implying that the 
time evolution of probability densities ceases to be tight, or the function N(t) may become a singular 
measure in finite time instead of being an Lj oc (R + ) function. Actually, in the numerical computations 
shown in Section (H we observe a blow-up in the value of the mean firing rate in finite time. This will 
need a modification of the notion of solution introduced in Definition 12.11 

Nevertheless, it is possible to obtain some a priori bounds with the help of appropriate choices of 
the test function <p m (|2.ip . Some of these choices are not allowed due to the growth at — oo of the 
test functions. We will say that a weak solution is fast-decaying at — oo if they are weak solutions 
in the sense of Definition 12.11 and the weak formulation in (|2.2p holds for all test functions growing 
algebraically in v. 

Lemma 2.3 (A priori estimates) Assume (jl.Tf) on the drift and diffusion coefficients and that 
(p,N) is a global-in-time solution of (|1.4p - ([1.6p in the sense of Definition \2.i\ fast decaying at —oo, 
then the following apriori estimates hold: 



(i) Ifb>V F -V R , then 





(ii) Ifb<VF-Vn then 




Moreover, if in addition a is constant then 
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Proof. With our decay assumption at — oo, we may use the test function (f>(v) = V F — v > 0. Then 
<[272|> gives 



d_ 

dt 

This is also written as 



V F pV F 

4>{v)p{v,t)dv = / [v-bN(t)]p{v,t)dv + N{t){V F -V R ). 



OO 



7 P Vf f 1 Vf 

- / tf>(v)p(v, t)dv+ / <f>(v)p(v, t)dv = V F - N(t) [b - (V F -Vr)}. (2.6) 

J — oo J — oo 

To prove (i), with our condition on b the term in N(t) is nonpositive and both results follow after 
integration in time. 

To prove (ii), we first use again (|2.6|) and, because the term in N(t) is nonnegative, we find the first 
result. Then, we use a truncation function <j)(v) £ C 2 such that 

<f>(V F ) = 1, <f>(v) = for v < V R , (j)'(v) > 0. 

Equation (12. 2D gives 
i rVp 

— 4>(v)p(v,t)dv + N(t) = 4>'(v)(- v + bN(t))p(v,t)dv + a <j)"(v)p(v, t) dv . 
dt J Vr J Vr ' J Vr 



Except regularity at v = Vr, we can have in mind (j)'(v) = 1/(V F — Vr) for v > Vr, and we can be as 
close as we want of this choice by paying a large second derivative of <j>. So with the parameter 



2 



1 



V f -Vr 



we may achieve with C(5) large 



— j (p(v)p(v, t) dv + 5N(t) <a <f>"(v)p(v, t) dv < C{5). 
dt J Vr J Vr 

And this leads directly to the result after integration in time. □ 

Corollary 2.4 Under the assumptions of Lemma [2.31 and assuming v 2 p°(v) E L x (— oo,Vf) and < 
b < V F — Vr, then the following apriori estimates hold: 

(i) If additionally a is constant, for all t > we have 

v 2 p(v,t)dv < C(l + t) 

oo 

(ii) If additionally -6min (v F , f^(V F - v)p°(v)dv)^ + o a + bV F + Yk^l < , th 

V F / r-V F 

v 2 p(v, t) dv < max I clq, / v 2 p°(v, t) dv 



en 
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Proof. We use <f>(v) = v 2 /2 as test function to get 



d f Vp v 2 



V F rVp y2 _ y2 

p(v,t)dv+ v 2 p(v,t)dv = bN(t) vp(v , t) dv + a(N (t)) + N (t) R F 



Vf 



lt, V 2 - V 2 

bV F + -2 £ 



V F \ V 2 - V 2 
6min ( V F , I (V F - v)p°(v)dv) )+a l + bV F + ^ 



= bN(t) / (v-V F )p(v,t)dv + a(N(t))+N(t) 

J —O 

<a + N(t) 

thanks to the first statement of Lemma 12.31 (ii). 

To prove (i), we just use the second statement of Lemma 12.31 (ii) valid for a constant which tells 
us that the time integration of the right-hand side grows at most linearly in time and so does 
v 2 p(v, t) dv. 

To prove (ii), we just use that the bracket is nonpositive and the results follows. □ 



3 Steady states 
3.1 Generalities 

This section is devoted to find all smooth stationary solutions of the problem (|1.4|) - f)1.6j) in the par- 
ticular relevant case of a drift of the form h(v) = Vq(N) — v. Let us search for continuous stationary 
solutions p of (jl.4p such that p is C 1 regular except possibly at V = Vr where it is Lipschitz. Using 
the definition in (12. 2D . we are then allowed by a direct integration by parts in the second derivative 
term of p to deduce that p satisfies 



d_ 

dv 



(v - V (N))p + a(N)—p(v) + NH(v - V R ) 
ov 







(3.1) 



in the sense of distributions, with H being the Heaviside function, i.e., H(u) = 1 for u > and 
H(u) = for u < 0. Therefore, we conclude that 

(v - V (N))p + a(N)^ + NH(v -V R ) = C. 
ov 



The definition of N in (jl.6p and the Dirichlet boundary condition f 1 1 . 5 j) imply C = by evaluating this 
expression at v = V F . Using again the boundary condition (]1.5p . p(V F ) = 0, we may finally integrate 
again and find that 



p(v) 



N (v-V (N)) 2 f V F (w-V (N)) 2 



a(N) 



e 2a H[w — Vji]dw 



which can be rewritten, using the expression of the Heaviside function, as 

,2 fVF 



p{v) 



N (v-V (N)) 

-. e I ( 

a \N ) Jmax(t),Vi{) 



(m-V Q (N))^ 



dw. 



(3.2) 
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Moreover, the firing rate in the stationary state N is determined by the normalization condition (|2.3[) 
or equivalently, 



a(N) 
N 



(v-V (N)Y 

e 2a 



V'f 



max(y,V R ) 



( W -v (N)y 



dw 



dv . 



(3.3) 



Summarizing all solutions with the above referred regularity of the stationary problem (|3,ip are of the 
form in (13. 2p with N being any positive solution to (j3. 3f) . 

Let us first comment that in the linear case Vq(N) = and a(N) = a > 0, we then get a unique 
stationary state given by the expression 



Poo{v) 



-e 2 ° 



e 2a dw. 



(3.4) 



max(u ,Vr) 



with iVoo the normalizing constant to unit mass over the interval (— oo, Vp], as obtained in [3]. 

The rest of this section is devoted to find conditions on the parameters of the model clarifying the 
number of solutions to (|3.3p . With this aim, it is convenient to perform a change of variables, and use 
new notations 



w-Vo 



Wf 



V F -Vo 



WR 



Vr-Vq 



(3.5) 



where the N dependency has been avoided to simplify notation. Then, we can rewrite the previous 
integral (and thus the condition for a steady state) as 

1 



N 



e 2 



tup 



e 2 du 



max(z,wn) 



(3.6) 



dz. 



Another alternative form of I(N) follows from the change of variables s = (z — u)/2 and s = (z + u)/2 
to get 

e' 2s ~ s ds ds = - 



I(N) = 2 
and consequently, 



e -2s* 



(e- 2sWF -e~ 2sWR ) ds. 



I(N) 



oo p -s 2 /2 



(e SWF -e SWR )ds. 



(3.7) 



3.2 Case of a(N) = a . 

We are now ready to state our main result on steady states. 

Theorem 3.1 Assume h(y, N) = b N — v , a(N) = a$ is constant and Vq = bN. 

i) For b < and b > small enough there is a unique steady state to (jl.4p - (jl.6p . 

ii) Under either the condition 



0<b<V F -V R , 



or the condition 



< 2a b < (V F - Vr) 2 Vr , 
then there exists at least one steady state solution to (|1.4j )- (|1.6p . 



(3.8) 
(3.9) 
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Figure 1: For several values of b, the function I(N) in (j3.6|) is plotted against the function 1/N (Left 
Figure) and the function NI(N) against the constant function 1 (Right Figure). Here a = 1, Vr = 1, 
V F = 2. 

Hi) If both (|3.9p and 6 > — Vr hold, then there are at least two steady states to (|1.4|) - (|1.6p . 
w) There is no steady state to (ll.4p - (|1.6p under the high connectivity condition 

b > max(2(V F - Vr), 2V f 1(0)). (3.10) 



Remark 3.2 It is natural to relate the absence of steady state for b large with blow-up of solutions. 
However, Theorem 12.21 in Section [2] shows this is not the only possible cause since the blow-up can 
happen for initial data concentrated enough around Vp independently of the value of b > 0. See also 
Section 5 for related numerical results. 

Proof. Let us first study properties of the function I(N). We first rewrite (|3.7p as 

I(N) = 



00 o/o -"eV^o - e v^o 
e, s ' z e. v^o 



ds. 



A direct Taylor expansion implies that 



e v^o — e v^o y F — Vr 



sVp 

< C s e 



(3.11) 



for all s > 0. Then, a direct application of the dominated convergence theorem and continuity theorems 
of integrals with respect to parameters show that the function I(N) is continuous on N on [0, oo). 
Moreover, the function I(N) is C°° on iV since all their derivatives can be computed by differentiating 
under the integral sign by direct application of dominated convergence theorems and differentiation 
theorems of integrals with respect to parameters. In particular, 



/'(AO 



oo Jo 



e - s Z/2, e sw F 



SW R 



)ds, 
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and for all integers k > 1, 

/( fe )(AT) 



■1) 



/a / Jo 
As a consequence, we deduce: 

1. Case b < 0: /(-/V) is an increasing strictly convex function and thus 

lim I(N) = oo . 

N^oo 

2. Case b > 0: /(iV) is a decreasing convex function. Also, it is obvious from the previous expansion 
(|3.11|) and dominated convergence theorem that 

lim I(N) = . 

N^oo 

It is also useful to keep in mind that, thanks to the form of I(N) in (|3.6|) . 

1(0) < v^K(O) - w R (0)]e max K(°),"!(0))/2 = V2^ {VF : VR) exp ( """Oft*® 
Now, let us show that for 6 > 0, we have 



"o 



2a 



< oo . (3.12) 



Using (]3.1ip . we deduce 



NI(N) - N 



V F -V R 



ao Jo 



lim N I(N) 

N^oo 



e -« / 2 e yso ^ s 



(3.13) 



< C iV 



se 1 e v^oev^ (is. 



A direct application of dominated convergence theorem shows that the right hand side converges to 
as ./V — > oo since sN exp(—^==) is a bounded function uniform in N and s. Thus, the computation of 
the limit is reduced to show 



lim N 

N^oo 



_ s 2 in son 

e v^o ds 



(3.14) 



With this aim, we rewrite the integral in terms of the complementary error function defined as 

2 f°° 2 

erfcix) := —= \ e~ l dt, 

V 71 " Jx 



and then 



00 _ s 2/ 2 _sWV b 2 N 2 fOO 

e v^o ds = e 2a o 



bN 



_( « + »£L)2 ^AF bijvi 

e V2 v 2a o ds = ^-=e 2a o erfc 

y/2 V\/2a 



Finally, we can obtain the limit (|3.14p using L'Hopital's rule 



lim N 



e 1 V^o ds 



TT 



V2 N 



lim 



erfc( 



bN ■ 
v / 2ao ' 



_ b 2 N 2 

e 2 a a 



V2 lim 



/2ao 



TV 



— — e 2a o 

a 



b 2 N 2 
.g 2a 
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With this analysis of the function I(N) we can now proof each of the statements of Theorem 13.11 



Proof of i). Let us start with the case b < 0. Here, the function I(N) is increasing, starting at 
1(0) < oo due to (|3.12p and such that 

lim I(N) = oo . 

Therefore, it crosses to the function 1/N at a single point. 

Now, for the case b > small, we first remark that similar dominated convergence arguments as 
above show that both I(N) and I'(N) are smooth functions of b. Moreover, it is simple to realize 
that I(N) is a decreasing function of the parameter b. Now, choosing < b < b* < (Vf — Vr)/2, then 
> I*(N) for all N > where I*(N) denotes the function associated to the parameter i>*. Using 
the limit (|3.13p . we can now infer the existence of iV* > depending only on b* such that 

NI(N) > NI*(N) > Vf ~~ Vr > 1 

for all N > N*. Therefore, by continuity of NI(N) there are solutions to NI(N) = 1 and all possible 
crossings of I(N) and 1/N are on the interval [0, iV*]. We observe that both I(N) and I'(N) converge 
towards the constant function 1(0) > and to respectively, uniformly in the interval [0, N*] as b — > 0. 
Therefore, for b small iV I(N) is strictly increasing on the interval [0, iV*] and there is a unique solution 
to N I(N) = 1. 



Proof of ii). Case of f|3.8jl . The claim that there are solutions to NI(N) = 1 for < b < Vp — Vr 
is a direct consequence of the continuity of I(N), A3. 12 j) and (j3.13j) . 

Case of ([3"U]) . We are going to prove that I(N) > 1/N for jy^^y2 < N < which concludes 
the existence of a steady state since 1(0) < oo due to f|3. 12|) implies that I(N) < 1/N for small N. 
Condition (|3.9jl only asserts that this interval for N is not empty. To do so, we show that 

(Vr - V F f 



I(N)> 



2rt 



for N G 



which obviously concludes the desired inequality I(N) > 1/N for the interval of iV under consideration. 

The condition > N is equivalent to wr > 0, therefore, using f|3 . 5f) and the expression for I(N) 
in (|3.6p . we deduce 



I(N)> 



e 2 du 



max(2,ji)fl) 



dz > 



W p 



VR 



e 2 



e 2 du 



dz . 



Since z > and e 2 is an increasing function for u > 0, then e 2 > e 2 on [z, wf], and we conclude 

(Vr - V F ? 



I(N)> 



Wp 



WR 



W F 



dudz 



2fto 



Proof of iii). Under the condition (|3,9p . we have shown in the previous point the existence of 
an interval where I(N) > I/N. On one hand, 1(0) < 00 in (|3.12|) implies that I(N) < I/N for N 
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small and the condition b > V F — Vr implies that I(N) < I/N for iV large enough due to the limit 
(|3.13j) . thus there are at least two crossings between I(N) and 1/JV. 

Proof of iv). Under assumption (|3.10p for b, it is easy to check that the following inequalities 
hold 

1(0) < 1/N for N < 2V F /b (3.15) 

and 

N>W F/i . (3.16) 

We consider N such that N > V F /b, this means that w F < 0. We use the formula (|3.7j) for I(N) and 
write the inequalities 



I(N) < (w F - w R ) / e~ s2/2 e SWF = (w F - w R )e w r /2 

Jo Jo 



e -(s-w F ) 2 /2 



>-w F \ W F\ y/a\w F \ 
where the mean-value theorem and w F < were used. Then, we conclude that 

/(AO < Vf ~ V * = Yl for N > V F /b. 
K ' sfa\w F \ bN — V F ' 

Therefore, using Inequality (j3. 16|) : 

I(N)<—, forN>2V F /b 

and due to the fact that / is decreasing and Inequality (j3. 15[) . we have I(N) < 1(0) < 1/N, for 
N < 2V F /b. In this way, we have shown that for all N, I(N) < 1/N and consequently there is no 
steady state. □ 



Remark 3.3 The functions I(N) and 1/N are depicted in Figured for the case Vq(N) = bN and 
a(N) = ao illustrating the main result: steady states exist for small b and do not exist for large b while 
there is an intermediate range of existence of two stationary states. The numerical plots of the function 
NI(N) might indicate that there are only three possibilities: one stationary state, two stationary states 
and no stationary state. However, we are not able to prove or disprove the uniqueness of a maximum 
for the function NI(N) eventually giving this sharp result. 

Remark 3.4 The condition (|3.9p can be improved by using one more term in the series expansion of 
the exponentials inside the integral of the expression of I(N) in (13.7f) . More precisely, ifw F > wr > 0, 
we use 

&SWF _ e s WR = ^i_ {w n_ w n ) ^Y^ S -(W F - W%). 
n=0 ' n=0 

In this way, we get 

T(m~> [°° -s*/2( Vf-Vr , If 2 2 ^ \ . ^ (V f -Vr)(V2^+(V f -Vr)) 
I(N)>J q e I \—— + -[ WF - WR ) s j ds> , 
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Figure 2: Left Figure: the function NI(N) against the constant 1 when a(N) is linear. For b = 0.5 we 
have considered a{N) = 0.5 + N/8, for b = 1.2: a(N) = 0.4 + iV/100 and for 6 = 8: a{N) = 6 + iV/100. 
Right Figure: the function NI(N) against the constant 1 with b < and a(N) = 1 + iV. Here Vr = 1; 
V F = 2. 



since 



/oo / poo 

e~ s2/2 ds = ^J-, J e~ s2/2 sds = l and V < Vr. 



Then, condition (|3.9p can be improved to 



2ab<V R {V F - V R ) + (V F - Vr) 



3.3 Case of a(N) = a + axN 

We now treat the case of a(N) = ao + aiN, with ao, ai > with 6 > 0. Proceeding as above we can 
obtain from (|3.7p the expression of its derivative 



d 



{h(N) - I 2 (N)) + — 



{V F h(N)-V R I 2 (N)), (3.17) 



where 



roc roc 

h(N)= / e- s2/2 e SWF ds and J 2 (iV) = / e~ s2 ' 2 e SWR ds . 
Jo Jo 

Therefore I(N) is decreasing since 



_d_ 

dN 



Vo(N) 
yfajN) 



2ba + ba 1 N 
2(a + ai iV) 3 / 2 



> and 



Ol 



dN \^/ajN) I (a + a!iV)3/2 



< 0. 
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Moreover, we can check that the computation of the limit ()3.13j) still holds. Actually, we have 



lim 



NI(N) 



N 



erfc 



bN 
\FTa 



T . T . lim ^ e -sV2 e - sb N/VE ds= lim ^ 



N 
v'Tn 



^erfc(ba) ^ 

/-V VrVl a — OOL . w . 



1 



a-s>oo -2fe 2 a 2 e-'' 2 a 2 _ e -h 2 a 2 5 ' 



where we have used the change a = and L'Hopital's rule. In the case b < 0, we can observe again 
by the same proof as before that I(N) — > oo when TV — > oo, and thus, by continuity there is at least 
one solution to NI(N) = 1. Nevertheless, it seems difficult to clarify perfectly the number of solutions 
due to the competing monotone functions in ([3.17]) . 

The generalization of part of Theorem 13. II is contained in the following result. We will skip its proof 
since it essentially follows the same steps as before with the new ingredients just mentioned. 

Corollary 3.5 Assume h(v, N) = b N — v, a(N) = ao + a\N with ao, oti > 0. 

i) Under either the condition b < Vf—Vr, or the conditions b > and 2ao&+2ai Vr < (Vf— Vr) 2 Vr, 
then there exists at least one steady state solution to (|1.4j) - (|1.6p . 

ii) If both 2a^b + 2a\VR < (Vf — Vr) Vr and b > Vf — Vr hold, then there are at least two steady 
states to ([L"4>([T1)]) . 

in) There is no steady state to ([Li]) - ([L6]) for b > max (2(V F - Vr),2V f 1(0)) . 

These behaviours are depicted in Figure [2j Let us point out that if a is linear and b < 0, I(N) may 
have a minimum for N > 0. 



4 Linear case and relaxation 

We study specifically the linear case, 6 = and a(N) = a, i.e., 

( dp^_ _ 9_ [vpM] _ ao ^p(v,t) = 5(v - V R )N(t), v < V F , 
p(V F ,t) = 0, N(t) := - ao -^p(V F , t) > 0, a > 0, 



p(v,Q) =p°(v) > 0, 



p°( V )dv = 1. 



(4.1) 



For later purposes, we remind that the steady state Poo(v) given in (I3.4|) satisfies 



d 2 

Poo(v) - 

Noo := -ao^ P oo(V F ) > 0, 



d 

-fol v Poo(v)] - a o-Q^2Poo(v) = 5(v - V R )N 00 , v < V F , 
Poo(Vf) = 0, 

V F 

Poo(v)dv = 1. 



(4.2) 
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We will assume in this section that solutions of the linear problem exist with the regularity needed in 
each result below and such that for all T > there exists Ct > such that p(v,t) < CtPoo{v) for all 
< t < T. These solutions might be obtained by the method developed in [10] and will be analysed 
elsewhere. 

We prove that the solutions converge in large times to the unique steady state Poo{v). Two relaxation 
processes are involved in this effect: dissipation by the diffusion term and dissipation by the firing 
term. This is stated in the following result about relative entropies for this problem. 

Theorem 4.1 Fast-decaying solutions to equation (|4,ip satisfy, for any smooth convex function G : 
R + — > R, the inequality 



d f Vp 

dt L P - (V)G 



f p(v,t) 

\Poo(v) 



,Poo{v) J V N °o Poo{v)j \Poo{v) 



." VF ( ^ n" ( P( V ' *) 



Poo{v) 



d ( p(v, t) 



dv \Pooiv) 



dv < 0. 



\v R 



(4.3) 



The following result is in fact standard on Poincare inequalities on R once q and p^, have been 
extended to the full line by odd symmetry with respect to Vf because p^ has a Gaussian behaviour 
at infinity thanks to (I3.4p . see [13]. 



Proposition 4.2 There exists v > such that 



f 

J — c 



V I Poo^V, 



< 



Vf 



Poo(v) 







q(v) 



dv \Poo(v) 



for all functions q such ^- G H 1 (p 00 (v)dv) . 

Note that performing the even symmetry of q with respect to Vf ensures that the extended function 
q satisfies 

q(v) dv = 0. 



These two theorems have direct consequences as in [15] . 

Corollary 4.3 (Exponential decay) Fast-decaying solutions to the equation (|4.1|) satisfy 



Vf 



, gMzgoM V < e -^ t [ Vf Poq{v) (p» -pM 



Poo(v) 



Poo(v) 



Proof. Taking q = p(v,t) — Poo(v) and G(x) = (x — l) 2 in the relative entropy inequality ()4.3p . we 
obtain 



d [ v * ,fp(v,t) 



2 r V F 
1 I < -2a 



/ Poo(v) 



d ( p(v, t) 



1 



dv \Poo(v) 

Poincare's inequality in Proposition 14.21 bounds the right hand side on the previous inequality 



a*i-oo \Poo{V) J J-oo \Poo(v) 

Finally, the Gronwall lemma directly gives the result. □ 
The proof of Theorem 14.11 is based on the following computations. 



1 
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Lemma 4.4 Given p a fast- decaying solution of (j4. 1|) . poo given by (|3.4p and £?(•) a convex function, 
then the following relations hold: 



d p 
dt p, 



v + A p ^) _ ao ^l^ = ^ S(v - V R ) (— - 

x \ Poo dv °°J dvpoo dv 2 Poo Poo KNoo Poo 



(4.4) 



Ot \Poo 



2ao d 

Poo dv 1 

P 

Poo 



dv^ ( 


' P 




,Poo 






dv Poo 





— ( 



p 

Poo 



Po 



N, 



00 Poo 



P 



(4.5) 



ub \l J oo 



d_ 

dv 



VPooG 



-OqPooG 



P 

Poo 

d p 

dv po 



a 



cP_ 

dv 2 



P00G 



P 



(4.6) 



-) + N eo 6(v-V R )\( 
Proof. Since l-f^W^I^-^-^we obtain 



Po 



Po 



dp d ( p \ p dp, 

TT- =Poo— I I + 

dv 



and 



d 2 p d 2 / p 

dv 2 dv 2 \p, 

Using these two expressions in 



dv VPoo / Poo dv 

d ( p \ dpoo , P d 2 p c 



Poo— I — I + 2— ( — ) + 



dv VPoo J dv poo dv 2 



(JL) = = J_ = v R )N(t) + %p{v,t)} + ao^p(M) 

VPoo / Poo at poo I <9?; 2 



we obtain 

Equation (|4.5p is a consequence of Equation (|4.4p and the following expressions for the partial 
derivatives of G ( — ) : 

( P_\ =G' ( — ( —\ —G ( — \ = G' ( — \ — (JL 

dt \Poo) \Poo) dt \Poo) ' dv \Poo) \Poo) dv \Poo 



and 



dv 2 \poo) VPoo/ \dv VPoo/ / VPoo/ dv 2 VPoo 



Finally, Equation (|4.6p is obtained using Equation (|4.5p and the fact that poo is solution of (|4.2p . □ 



Proof Theorem 14.11 We integrate from —00 to Vf — a in (|4.6|) and let a tend to + and use 
L'Hopital's rule 

p(v,t) 



lim 



lim 



(4.7) 
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Since p(v, t) < CtPoo with < t < T, then 



d_ 

dt 



p \ d 
PooG ( — ) dv - a — 

K Poo ) OV 



PooG 



-ao 



V F 



PooG"! — ] (^— ) dv + N t 



Po 



d p 



dv p c 



N p 



Noo Po 



P 



Vr 



The Dirichlet boundary condition (|1.5|) implies that 



d_ 

dv 



PooG 



p_ 

Poc 



\V F 



-ao- 



dpo 



dv 



where we used that 

9 n ( P 
Poo-^-G — 

dv \p a 



PooG' 



, ( p \ -Npoo + A^oo p 



Po 



Plo a o 



v v = G 



q 

Poo 



-N.N^p 

1 Vf - U ' 

«0 «0 Poo 



due to (|4.7p . Collecting all terms leads to the desired inequality. 



□ 



5 Numerical results 

We consider an explicit method to simulate the numerical approximation for the NNLIF ()1.4|) . We 
base our algorithm on standard shock-capturing methods for the advection term and second-order 
finite differences for the second-order term. More precisely, the first order term is approximated by 
finite difference WENO-schemes |21| . 




Figure 3: Distribution functions p(v, t) for b = 0.5 and a = 1 at differents times. 



The time evolution is performed with a TVD Runge-Kutta scheme. Other finite difference scheme 
for the Fokker-Planck equation has been used such as the Chang-Cooper method [5j similarly to the 
approximation studied in [6] for a model with variable voltage and conductance. The Chang-Cooper 
method presents difficulties when the firing rate becomes large and the diffusion coefficient a(N) is 
constant. To discuss this, we have just to remind the reader that the Chang-Cooper method performs 
a kind of ^-finite difference approximation of p/M where M is a Maxwellian in the kernel of the 
linear Fokker-Planck operator. Whenever a(N) is constant, b > and N is large, the drift of the 
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Maxwellian, in terms of which is rewritten the Fokker-Planck equation, practically vanishes on the 
interval (— oo, Vf] and this particular Chang-Cooper method is not suitable. 

In our simulations we consider a uniform mesh in v, for v £ {Ymim^F]- The value V m in (less than 
Vr) is adjusted in the numerical experiments to fulfill that p(V m i n ,t) ~ 0, while Vf is fixed to 2 and 
Vr = 1. Most of our initial data are Maxwellians: 

^ (d-Dq) 2 

Po(v) = j=e 2CT o , 

0"oV 27T 

where the mean vq and the variance Oq are chosen according to the analyzed phenomenon. When 
the system has two steady states, we also take as initial data the profiles given by ()3.2p with N an 
approximate value of the stationary firing rate, in order to start close to the stationary state with 
larger firing rate. 




Figure 4: Firing rates N(t) for a = 1. Top left: b = 0.5 with initial data a Maxwellian with: vq = 
and o"q = 0.25. Top right: 6 = 3 with initial data a Maxwellian with: vo = —1 and = 0.5. Bottom 
left: b = 1.5 considering two different initial data: a Maxwellian with: vo = —1 and = 0.5 and a 
profile given by the expression (|3,2p with iV = 2.31901. Bottom right: b = —1.5 with initial data a 
Maxwellian with: vo = —1 and = 0.5. The top right case seems to depict a blow-up phenomena 
demonstrated in Theorem 12.21 

Steady states.- As we show in Section [31 for b positive there is a range of values for which there are 
either one or two or no steady states. With our simulations we can observe all the cases represented 
in Figures [Q and [2j 
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Figure 5: For b = 1.5 and a(N) = 1 figures show unstability of the steady state with higher firing rate. 
Left: Evolution on time of the firing rate considering different initial firing rate. Right: Evolution on 
time of the distribution function with initial firing rate 2.31901. In both figures we have considered 
V R = 1; V F = 2. 



In Figure Owe show the time evolution of the distribution function p(v, t), in the case of a = 1 and 
b = 0.5, considering as initial data a Maxwellian with vq = and <Tq = 0.25. We observe that the 
solution after 3.5 time units numerically achieves the steady state with the imposed tolerance. The top 
left subplot in Figure 2] describes the time evolution of the firing rate, which becomes constant after 
some time. This clearly corresponds to the case of a unique locally asymptotically stable stationary 
state. Let us remark that in the right subplot of Figure [3j we can observe the Lipschitz behavior of 
the function at Vr as it should be from the jump in the flux and thus on the derivative of the solutions 
and the stationary states, see Section [3j 

For b = 1.5, we proved in Section [3] that there are two steady states. With our simulations we can 
conjecture that the steady state with larger firing rate is unstable. However the stationary solution 
with low firing rate is locally asymptotically stable. We illustrate this situation in the bottom left 
subplot in Figure [H Starting with a firing rate close to the high stationary firing value, the solution 
tends to the low firing stationary value. 

In Figure [5] we analyze in more details the behavior of the steady state with larger firing rate. The 
left subplot presents the evolution on time of the firing rate for different distribution function starting 
with profiles given by the expression (|3.2j) with iV an approximate value of the stationary firing rate. 
We show that, depending of the initial firing rate considered, its behavior is different: tends to the 
lower steady state or goes to infinity. The firing rate for the solution with initial Nq = 2.31901 remains 
almost constant for a period of time. Observe in Figure [5] that the difference between the initial data 
and the distribution function at time t = 1.8 is almost negligible. However, the system evolves slowly 
and at t = 6 the distribution is very close to the the lower steady state, see the bottom left subplot in 
Figure [U 

In the bottom right subplot of Figure H] we observe the evolution for a negative value of b, where we 
know that there is always a unique steady state, and its local asymptotic stability seems clear from 
the numerical experiments. 
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No steady states.- Our results in Section [3] indicate that there are no steady states for b = 3. In 
Figure [6] we observe the evolution on time of the distribution function p. In Figure U] (right top) we 
show the time evolution of the firing rate, which seems to blow up in finite time. We observe how the 
distribution function becomes more and more picked at Vr and Vf producing an increasing value of 
the firing rate. 

Blow up.- According to our blow-up Theorem 12.21 the blow-up in finite time of the solution happens 
for any value of b > if the initial data is concentrated enough on the firing rate. In Figures [7] and [HI 
we show the evolution on time of the firing rate with an initial data with mass concentrated close to 
Vf for values of b in which there are either a unique or two stationary states. The firing rate increases 
without bound up to the computing time. It seems that the blow-up condition in Theorem 12.21 is not 
as restrictive as to say that the initial data is close to a Dirac Delta at Vf- Let us finally mention 
that blow-up appears numerically also in case of a(N) = a^ + a±N, but here the blow-up scenario is 
characterized by a break-up of the condition under which (jl.6p has a unique solution N, i.e., 



§>■<> 



< 1 



Therefore, the blow-up in the value of the firing rate appears even if the derivative of p at the firing 
voltage does not diverge. 



6 Conclusion 

The nonlinear noisy leaky integrate and fire (NNLIF) model is a standard Fokker-Planck equation 
describing spiking events in neuron networks. It was observed numerically in various places, but 
never stated as such, that a blow-up phenomena can occur in finite time. We have described a class 
of situations where we can prove that this happens. Remarkably, the system can blow-up for all 
connectivity parameter b > 0, whatever is the (stabilizing) noise. 

The nature of this blow-up is not mathematically proved. Nevertheless, our estimates in Lemma 
12.31 indicate that it should not come from a vanishing behaviour for !) ~ -oo, or a lack of fast decay 
rate because the second moment in v is controlled uniformly in blow-up situations. Additionally, 
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Figure 7: Parameter values are a = 1 and 6 = 1.5 and this corresponds to two steady states. Left: 
Evolution of the distribution function p(v,t) in time of an initial Maxwellian centered at v = 1.5 and 
with variance 0.005. Right: Time evolution of the firing rate; again we observe numerically a blow-up 
behaviour for an initial data enough concentrated near Vp. 

numerical evidence is that the firing rate N(t) blows-up in finite time whenever a singularity in the 
system occurs. This scenario is compatible with all our theoretical knowledge on the NNLIF and in 
particular with L 1 estimates on the total network activity (firing rate N(t)). 

Blow-up has also been proved to occur in the deterministic quadratic adaptive Integrate-and-Fire 
model in |22j . The blow-up scenario here is quite different from ours in many aspects; the model is a 
linear model (in our terminology) for a single neuron not a network. Remarkably, the blow-up scenario 
arises on the adaptation variable that we do not have in the NNLIF. These are interesting questions 
to know if blow-up could occur in a noisy 'linear' adaptive situation. 

We have established that the set of steady states can be empty, a single state or two states depending 
on the network connectivity. These are all compatible with blow-up profile, and when they exist, 
numerics can exhibit convergence. Several questions are left open; is it possible to have triple or more 
steady states? Which of them are stable? Can a bifurcation analysis help to understand and compute 
the set of steady states? 
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Evolution of the distribution function p(v, t) in time of an initial Maxwellian centered at v = 1.83 and 
with variance 0.003. Right: Time evolution of the firing rate; again this seems to be a typical blow-up 
behaviour. 
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